package com.salinity.kun.algorithm.harmonic;

import java.util.ArrayList;
import java.util.Arrays;
import java.util.List;
import java.util.Map;
import com.salinity.kun.helper.BaseHydrologyDataHelper;
import com.salinity.kun.helper.MatrixHelper;
import com.salinity.kun.model.StationDateData;

import Jama.Matrix;

public class Harmonic {

	private List<StationDateData> dataList;
	private double[] tides;
	private double deltaT;// 小时
	private List<TideSymbol> tideSymbolList;
	private Map<TideSymbol, Double> angelMap;

	/**
	 * 
	 * @param dataList
	 * @param deltaT
	 *            单位为小时
	 */
	public Harmonic(List<StationDateData> dataList, double deltaT) {
		this.dataList = dataList;
		this.tides = BaseHydrologyDataHelper.hydToDoubleList(dataList);
		this.deltaT = deltaT;
	}

	private Map<TideSymbol, Double> getAngelMap() {
		return angelMap == null ? initAngelMap() : angelMap;
	}

	private double getAngel(TideSymbol tideSymbol) {
		return getAngelMap().get(tideSymbol) * Math.PI / 180;
	}

	private List<TideSymbol> getCalTideSymbolList() {
		return tideSymbolList == null ? initTideList() : tideSymbolList;
	}

	public void setTideSymbolList(List<TideSymbol> tideSymbolList) {
		this.tideSymbolList = tideSymbolList;
	}

	private Map<TideSymbol, Double> initAngelMap() {
		angelMap = ConstantParamHelper.getAngularVelocity().getAngel();
		return angelMap;
	}

	private List<TideSymbol> initTideList() {
		tideSymbolList = new ArrayList<>();
		tideSymbolList.add(TideSymbol.M2);
		tideSymbolList.add(TideSymbol.S2);
	    tideSymbolList.add(TideSymbol.N2);
	    tideSymbolList.add(TideSymbol.K2);

		tideSymbolList.add(TideSymbol.K1);
		tideSymbolList.add(TideSymbol.O1);
		 tideSymbolList.add(TideSymbol.P1);
		tideSymbolList.add(TideSymbol.Q1);

		tideSymbolList.add(TideSymbol.M4);	
	

		tideSymbolList.add(TideSymbol.M6);
		tideSymbolList.add(TideSymbol.MS4);
		// tideSymbolList.add(TideSymbol.M2);
		// tideSymbolList.add(TideSymbol.S2);
		// tideSymbolList.add(TideSymbol.N2);
		// tideSymbolList.add(TideSymbol.K1);
		// tideSymbolList.add(TideSymbol.O1);
		// tideSymbolList.add(TideSymbol.Q1);
		// tideSymbolList.add(TideSymbol.M4);
		// tideSymbolList.add(TideSymbol.MS4);
		// tideSymbolList.add(TideSymbol.M6);

		return tideSymbolList;
	}

	/**
	 * 
	 * replaced by 'checkRslt(double[] ax, double[] bx, double ti)'
	 * 
	 * @param p
	 * @param t
	 *            单位为：hour
	 * @return
	 */
	@Deprecated
	public double checkRslt(double[] p, double ti) {
		List<TideSymbol> tideSymbolList = getCalTideSymbolList();
		double rslt = 0, W;
		int n = (p.length - 1) / 2;
		rslt += p[0];

		for (int i = 1; i <= n; i++) {
			W = getAngel(tideSymbolList.get(i - 1));
			rslt += p[i] * Math.cos(W * ti) + p[i + n] * Math.sin(W * ti);
		}
		return rslt;
	}

	/**
	 * 根据参数和时间t计算 所得值，用于与原值进行比较
	 * 
	 * @param ax
	 * @param bx
	 * @param ti
	 *            单位为：hour
	 * @return
	 */
	public double checkRslt(double[] ax, double[] bx, double ti) {
		double rslt = 0, W;
		List<TideSymbol> tideSymbolList = getCalTideSymbolList();
		rslt += ax[0];
		int i;
		for (i = 1; i < bx.length; i++) {
			W = getAngel(tideSymbolList.get(i - 1));
			rslt += ax[i] * Math.cos(W * ti) + bx[i - 1] * Math.sin(W * ti);
		}
		return rslt;
	}

	/**
	 * 周期为[-T/2,T/2]
	 * 
	 *
	 * @param offset_t
	 * 
	 *            起点时刻与最开始起点的偏移小时数（如当前时间起点与2017年1月1日 0：0:0的偏移，正数）
	 * 
	 * @return 结果分别为:a,fa, 第二个为b，fb
	 */
	public List<double[][]> getMatrix_v1(double offset_t) {
		List<TideSymbol> tideSymbolList = getCalTideSymbolList();
		double ti, w, ds1, ds2;
		int i, j, k;
		int m = tideSymbolList.size();

		int len = tides.length;
		double L = (len - 1) * deltaT;
		double[][] a = new double[m + 1][m + 1];
		double[][] fa = new double[m + 1][1];
		double[][] b = new double[m][m];
		double[][] fb = new double[m][1];

		// [-L/2,L/2]
		double halfT = L / 2f;
		double real_t0 = offset_t - halfT;
		double real_tn = offset_t + halfT;

		// a[0][0]
		a[0][0] = L;

		// 求 fa(0),这里是通过求梯形面积和来实现对（nt） 积分
		fa[0][0] = 0.5 * (tides[0] + tides[len - 1]) * deltaT;
		for (i = 1; i < len - 1; i++) {
			fa[0][0] += tides[i] * deltaT;
		}

		// 求 fa(1~m)
		for (i = 1; i < m + 1; i++) {
			w = getAngel(getCalTideSymbolList().get(i - 1));

			// fa(1), 始末两项分别为 n(0)cos(-L/2) 、n(len-1)cos(L/2),cos为偶函数，故 cos(-L/2) = cos(L/2)
			// fa[i][0] = 0.5 * (tides[0] + tides[len - 1]) * Math.cos(w * t0) * deltaT;
			fa[i][0] = 0.5 * (tides[0] * Math.cos(w * real_t0) + tides[len - 1] * Math.cos(w * real_tn)) * deltaT;
			// fb(1)
			// fb[i - 1][0] = 0.5 * (-tides[0] + tides[len - 1]) * Math.sin(w * t0) *
			// deltaT;
			fb[i - 1][0] = 0.5 * (tides[0] * Math.sin(w * real_t0) + tides[len - 1] * Math.sin(w * real_tn)) * deltaT;
			for (k = 1; k < len - 1; k++) {
				ti = real_t0 + deltaT * k;

				// fa
				fa[i][0] += tides[k] * Math.cos(w * ti) * deltaT;

				// fb
				fb[i - 1][0] += tides[k] * Math.sin(w * ti) * deltaT;
			}
		}

		// 求ax,bx
		for (i = 1; i < m + 1; i++) {
			w = getAngel(getCalTideSymbolList().get(i - 1));
			a[i][0] = 2.0 * Math.sin(w * halfT) / w;
			for (j = 1; j < m + 1; j++) {
				ds1 = getAngel(getCalTideSymbolList().get(j - 1)) - getAngel(getCalTideSymbolList().get(i - 1));
				ds2 = getAngel(getCalTideSymbolList().get(j - 1)) + getAngel(getCalTideSymbolList().get(i - 1));
				if (j != i) {

					// fa
					a[i][j] = Math.sin(ds2 * halfT) / ds2 + Math.sin(ds1 * halfT) / ds1;

					// fb
					b[i - 1][j - 1] = Math.sin(ds1 * halfT) / ds1 - Math.sin(ds2 * halfT) / ds2;

				} else {

					// fa, 此时ds2 = 2sigma(j),则sin(sigma(j)*L)=sin(2sigma(j)*(L/2))=sin(ds2*t0)
					a[i][j] = halfT + Math.sin(ds2 * halfT) / ds2;

					// fb
					b[i - 1][j - 1] = halfT - Math.sin(ds2 * halfT) / ds2;
				}
			}
		}

		// 补全第一行，因为对称，直接将第一列赋值给第一行
		for (k = 1; k <= m; k++) {
			a[0][k] = a[k][0];
		}

		List<double[][]> rslt = new ArrayList<>();
		rslt.add(a);
		rslt.add(fa);
		rslt.add(b);
		rslt.add(fb);
		return rslt;
	}

	/**
	 * 周期为[-T/2,T/2]
	 * 
	 * 
	 * @return 结果分别为:a,fa, 第二个为b，fb
	 */
	public List<double[][]> getMatrix_v1() {
		List<TideSymbol> tideSymbolList = getCalTideSymbolList();
		double ti, w, ds1, ds2;
		int i, j, k;
		int m = tideSymbolList.size();

		int len = tides.length;
		double L = (len - 1) * deltaT;
		double[][] a = new double[m + 1][m + 1];
		double[][] fa = new double[m + 1][1];
		double[][] b = new double[m][m];
		double[][] fb = new double[m][1];

		// [-L/2,L/2]
		double halfT = L / 2f;
		double real_t0 = -halfT;
		double real_tn = halfT;

		// a[0][0]
		a[0][0] = L;

		// 求 fa(0),这里是通过求梯形面积和来实现对（nt） 积分
		fa[0][0] = 0.5 * (tides[0] + tides[len - 1]) * deltaT;
		for (i = 1; i < len - 1; i++) {
			fa[0][0] += tides[i] * deltaT;
		}

		// 求 fa(1~m)
		for (i = 1; i < m + 1; i++) {
			w = getAngel(getCalTideSymbolList().get(i - 1));

			// fa(1), 始末两项分别为 n(0)cos(-L/2) 、n(len-1)cos(L/2),cos为偶函数，故 cos(-L/2) = cos(L/2)
			// fa[i][0] = 0.5 * (tides[0] + tides[len - 1]) * Math.cos(w * t0) * deltaT;
			fa[i][0] = 0.5 * (tides[0] * Math.cos(w * real_t0) + tides[len - 1] * Math.cos(w * real_tn)) * deltaT;
			// fb(1)
			// fb[i - 1][0] = 0.5 * (-tides[0] + tides[len - 1]) * Math.sin(w * t0) *
			// deltaT;
			fb[i - 1][0] = 0.5 * (tides[0] * Math.sin(w * real_t0) + tides[len - 1] * Math.sin(w * real_tn)) * deltaT;
			for (k = 1; k < len - 1; k++) {
				ti = real_t0 + deltaT * k;

				// fa
				fa[i][0] += tides[k] * Math.cos(w * ti) * deltaT;

				// fb
				fb[i - 1][0] += tides[k] * Math.sin(w * ti) * deltaT;
			}
		}

		// 求ax,bx
		for (i = 1; i < m + 1; i++) {
			w = getAngel(getCalTideSymbolList().get(i - 1));
			a[i][0] = 2.0 * Math.sin(w * halfT) / w;
			for (j = 1; j < m + 1; j++) {
				ds1 = getAngel(getCalTideSymbolList().get(j - 1)) - getAngel(getCalTideSymbolList().get(i - 1));
				ds2 = getAngel(getCalTideSymbolList().get(j - 1)) + getAngel(getCalTideSymbolList().get(i - 1));
				if (j != i) {

					// fa
					a[i][j] = Math.sin(ds2 * halfT) / ds2 + Math.sin(ds1 * halfT) / ds1;

					// fb
					b[i - 1][j - 1] = Math.sin(ds1 * halfT) / ds1 - Math.sin(ds2 * halfT) / ds2;

				} else {

					// fa, 此时ds2 = 2sigma(j),则sin(sigma(j)*L)=sin(2sigma(j)*(L/2))=sin(ds2*t0)
					a[i][j] = halfT + Math.sin(ds2 * halfT) / ds2;

					// fb
					b[i - 1][j - 1] = halfT - Math.sin(ds2 * halfT) / ds2;
				}
			}
		}

		// 补全第一行，因为对称，直接将第一列赋值给第一行
		for (k = 1; k <= m; k++) {
			a[0][k] = a[k][0];
		}

		List<double[][]> rslt = new ArrayList<>();
		rslt.add(a);
		rslt.add(fa);
		rslt.add(b);
		rslt.add(fb);
		return rslt;
	}

	/**
	 * 周期为[0,T]
	 * 
	 * @return 结果为ab,f(ab)
	 */
	public List<double[][]> getMatrix_v2() {
		List<double[][]> rslt = new ArrayList<>();
		List<TideSymbol> tideSymbolList = getCalTideSymbolList();
		int m = tideSymbolList.size();
		double[][] fab = new double[2 * m + 1][1];
		double[][] ab = new double[2 * m + 1][2 * m + 1];

		int i, j, k;
		double w, ti, ds1, ds2;
		int len = tides.length;
		double T = (len - 1) * deltaT;

		// ab[0][0]
		ab[0][0] = T;

		// 计算 fab(0), 通过求梯形面积和来实现对（nt） 积分
		fab[0][0] = 0.5 * (tides[0] + tides[len - 1]) * deltaT;

		for (i = 1; i < len - 1; i++) {
			fab[0][0] += tides[i] * deltaT;
		}

		// 求 fab(1~m)
		for (i = 1; i < m + 1; i++) {
			w = getAngel(getCalTideSymbolList().get(i - 1));

			// 始末两项分别为 n(0)cos(w*0) 、n(len-1)cos(W*T);
			fab[i][0] = 0.5 * (tides[0] + tides[len - 1] * Math.cos(w * T)) * deltaT;

			// 始末两项分别为n(0)sin(w*0) 、n(len-1)sin(W*T);
			fab[i + m][0] = 0.5 * (tides[len - 1] * Math.sin(w * T)) * deltaT;

			for (k = 1; k < len - 1; k++) {
				ti = 0.0 + deltaT * k;

				fab[i][0] += tides[k] * Math.cos(w * ti) * deltaT;

				fab[i + m][0] += tides[k] * Math.sin(w * ti) * deltaT;
			}
		}

		// 求 ab
		for (i = 1; i < m + 1; i++) {
			w = getAngel(getCalTideSymbolList().get(i - 1));

			ab[i][0] = Math.sin(w * T) / w;
			ab[i + m][0] = (1 - Math.cos(w * T)) / w;

			for (j = 1; j < m + 1; j++) {
				ds1 = getAngel(getCalTideSymbolList().get(j - 1)) - getAngel(getCalTideSymbolList().get(i - 1));
				ds2 = getAngel(getCalTideSymbolList().get(j - 1)) + getAngel(getCalTideSymbolList().get(i - 1));

				if (i != j) {

					// a(i)(1~m)的系数
					ab[i][j] = 0.5 * (Math.sin(ds1 * T) / ds1 + Math.sin(ds2 * T) / ds2);

					// b(i)(1~m)的系数
					ab[i][j + m] = 0.5 * ((1 - Math.cos(ds1 * T)) / ds1 + (1 - Math.cos(ds2 * T)) / ds2);

					// a(i+m)(1~m)的系数,注意ds1前面的负号！！！！
					ab[i + m][j] = 0.5 * ((1 - Math.cos(-ds1 * T)) / -ds1 + (1 - Math.cos(ds2 * T)) / ds2);

					// b(i+m)(1~m)的系数
					ab[i + m][j + m] = 0.5 * (Math.sin(ds1 * T) / ds1 - Math.sin(ds2 * T) / ds2);

				} else {

					// a(i)(1~m)的系数
					ab[i][j] = 0.5 * (T + Math.sin(2 * w * T) / (2 * w));

					// b(i)(1~m)的系数
					ab[i][j + m] = Math.sin(w * T) * Math.sin(w * T) / (2 * w);

					// a(i+m)(1~m)的系数,因为对称，直接赋值
					ab[i + m][j] = Math.sin(w * T) * Math.sin(w * T) / (2 * w);

					// b(i+m)(1~m)的系数
					ab[i + m][j + m] = 0.5 * (T - Math.sin(2 * w * T) / (2 * w));
				}
			}

		}

		// 补全第一行，因为对称，直接将第一列赋值给第一行
		for (k = 1; k < 2 * m + 1; k++) {
			ab[0][k] = ab[k][0];
		}
		handleParamsMatrix(ab);
		rslt.add(ab);

		rslt.add(fab);

		return rslt;
	}

	/**
	 * 对系数矩阵中共线的项进行处理
	 */
	public double[][] handleParamsMatrix(double[][] m) {
		// 这里直接处理 K2(2,4)，P1(5,7) ,{ 2 + rows, 4 + cols }, { 5 + rows, 7 + cols }

		int rows = m.length / 2;
		int cols = m[0].length / 2;

		double[][] rslt = new Matrix(m).getArray();

		int handleIndex[][] = { { 2, 4 }, { 5, 7 }, { 2 + rows, 4 + cols }, { 5 + rows, 7 + cols } };

		int i, j, k;
		// MatrixHelper.printMatrix(rslt);
		for (k = 0; k < handleIndex.length; k++) {
			i = handleIndex[k][0];
			j = handleIndex[k][1];
			// rslt[i][j]*=0.996;
			// rslt[j][i]*=0.996;
		}

		// MatrixHelper.printMatrix(rslt);
		return rslt;

	}

	/**
	 * 
	 * @param m
	 *            系数矩阵 ax=b 中的a
	 * @param b
	 *            ax=b 中的b
	 * @param offset
	 *            岭回归偏移系数
	 * @return 第一个为ax的解,第二个为bx的解
	 */
	public static Matrix calParams(double[][] m, double[][] b, double offset) {
		Matrix A = new Matrix(m);
		Matrix Ab = new Matrix(b);

		Ab = Ab.times(1 / A.get(0, 0));
		A = A.times(1 / A.get(0, 0));

		// 求ax x向量 = c * () * c^T*b
		Matrix lamda = A.eig().getD();

		 //System.out.println(A.cond());
		
//		 MatrixHelper.printMatrix(A.getArray(), "\t");
//		 
//		MatrixHelper.printMatrix(Ab.transpose().getArray(), "\t");
//
//		
//		MatrixHelper.printMatrix(lamda.getArray(), "\t");
//		System.out.println();
		Matrix im = Matrix.identity(lamda.getRowDimension(), lamda.getColumnDimension());

		// 利用岭回归的思想，调整lamda
		double step = offset;
		Matrix c = A.eig().getV();
		
		c.times(c.transpose()).print(4,6);
		
		Matrix ic = c.inverse();
		Matrix iLamda, X;

		iLamda = lamda.plus(im.times(step)).inverse();
		// MatrixHelper.printMatrix(iLamda.getArray(),",");
		X = c.times(iLamda).times(ic).times(Ab);
		return X;

	}

	/**
	 * 
	 * @param m
	 * @param b
	 * @param threshold
	 *            只对小于阈值的数进行变化
	 * @param offset
	 * @return
	 */
	public static Matrix calParams(double[][] m, double[][] b, double threshold, double offset) {
		Matrix A = new Matrix(m);
		Matrix Ab = new Matrix(b);

		Ab = Ab.times(1 / A.get(0, 0));
		A = A.times(1 / A.get(0, 0));

		// 求ax x向量 = c * () * c^T*b
		Matrix lamda = A.eig().getD();

		Matrix lamda1 = handleLingParams(lamda, threshold, offset);

		// lamda1.print(4,6);

		// 利用岭回归的思想，调整lamda
		double step = offset;
		Matrix c = A.eig().getV();
		Matrix ic = c.inverse();
		Matrix iLamda, X;

		iLamda = lamda.plus(lamda1.times(step)).inverse();

		X = c.times(iLamda).times(ic).times(Ab);
		return X;

	}

	public static Matrix handleLingParams(Matrix lamda, double threshold, double offset) {
		int i, m = lamda.getRowDimension(), n = lamda.getColumnDimension();
		Matrix rslt = new Matrix(lamda.getArray());
		for (i = 0; i < m; i++) {
			if (lamda.get(i, i) < threshold) {
				rslt.set(i, i, lamda.get(i, i) + offset);
			}
		}
		return rslt;
	}

	public static Matrix calParams(double[][] matrix, double[][] b) {
		return calParams(matrix, b, 0.001);
	}

	/**
	 * 将对角元素化为1
	 */
	public static double[][] normalizeMatrix(double[][] m) {
		int i, j, row = m.length, col = m[0].length;
		double temp;
		double[][] rslt = new double[row][col];

		for (i = 0; i < row; i++) {
			temp = m[i][i];
			for (j = 0; j < col; j++) {
				rslt[i][j] = m[i][j] / temp;
			}
		}
		return rslt;
	}

	public static double[][] matrixDivNum(double[][] m, double num) {
		int i, j, row = m.length, col = m[0].length;
		double[][] rslt = new double[row][col];

		for (i = 0; i < row; i++) {
			for (j = 0; j < col; j++) {
				rslt[i][j] = m[i][j] / num;
			}
		}
		return rslt;
	}

	/**
	 * 用此方法容易产生病态矩阵 replace by getMatrix_v1
	 * 
	 * @return
	 */
	@Deprecated
	public double[][] getMatrix() {
		List<TideSymbol> tides = getCalTideSymbolList();
		int rows = tides.size() * 2 + 1;
		int cols = rows + 1;
		int n = dataList.size();
		double[][] rslt = new double[rows][cols];
		double hi;
		double ti = 0.0;
		for (int k = 0; k < n; k++) {
			hi = dataList.get(k).getValue();
			ti = deltaT * k;

			calParams(rslt, tides, hi, ti);
			// GaussianElimination.printMatrix(rslt);
			// System.out.println("----------");
		}
		return rslt;
	}

	public static void divideT(double[][] m, double T) {
		int rows = m.length;
		int cols = m[0].length;

		for (int i = 0; i < rows; i++) {
			for (int j = 0; j < cols; j++) {
				m[i][j] /= T;
			}
		}
	}

	@Deprecated
	public void calParams(double[][] rslt, List<TideSymbol> tides, double hi, double ti) {
		int i, j, J = tides.size(), rows = 2 * J + 1, cols = rows + 1;
		double Wi, Wj;
		// 最左上角,每次加1
		rslt[0][0] += 1;

		// 最右上角，hi
		rslt[0][cols - 1] += hi;

		// 第一列上半部分：cos(Wj*ti)
		for (i = 1; i <= J; i++) {
			Wi = getAngel(tides.get(i - 1));

			rslt[i][0] += Math.cos(Wi * ti);
		}

		// 第一列下半部分：sin(Wj*ti)
		for (; i <= 2 * J; i++) {
			Wi = getAngel(tides.get(i - 1 - J));
			rslt[i][0] += Math.sin(Wi * ti);
		}

		// 最后一列,上半部分 cos(Wj*ti)*hi
		for (i = 1; i <= J; i++) {
			Wi = getAngel(tides.get(i - 1));
			rslt[i][cols - 1] += Math.cos(Wi * ti) * hi;
		}

		// 最后一列,下半部分 sin(Wj*ti)*hi
		for (; i <= 2 * J; i++) {
			Wi = getAngel(tides.get(i - 1 - J));
			rslt[i][cols - 1] += Math.sin(Wi * ti) * hi;
		}

		// 第一行除开开头结尾,左半部分：cos(Wj*ti)
		for (j = 1; j <= J; j++) {
			Wj = getAngel(tides.get(j - 1));
			rslt[0][j] += Math.cos(Wj * ti);
		}

		// 第一行除开开头结尾,右半部分：sin(Wj*ti)
		for (; j <= 2 * J; j++) {
			Wj = getAngel(tides.get(j - 1 - J));
			rslt[0][j] += Math.sin(Wj * ti);
		}
		// 中间左上 cos(Wj*ti)cos(Wi*ti)
		for (i = 1; i <= J; i++) {
			Wi = getAngel(tides.get(i - 1));
			for (j = 1; j <= J; j++) {
				Wj = getAngel(tides.get(j - 1));
				rslt[i][j] += Math.cos(Wj * ti) * Math.cos(Wi * ti);
			}
		}

		// 中间左下cos(Wj*ti)sin(Wi*ti)
		for (; i <= 2 * J; i++) {
			Wi = getAngel(tides.get(i - 1 - J));
			for (j = 1; j <= J; j++) {
				Wj = getAngel(tides.get(j - 1));
				rslt[i][j] += Math.cos(Wj * ti) * Math.sin(Wi * ti);
			}
		}

		// 中间右上 sin(Wj*ti)*cos(Wi*ti)
		for (i = 1; i <= J; i++) {
			Wi = getAngel(tides.get(i - 1));
			for (j = J + 1; j <= 2 * J; j++) {
				Wj = getAngel(tides.get(j - 1 - J));
				rslt[i][j] += Math.sin(Wj * ti) * Math.cos(Wi * ti);
			}
		}

		// 中间右下sin(Wj*ti)*sin(Wi*ti)
		for (; i <= 2 * J; i++) {
			Wi = getAngel(tides.get(i - 1 - J));
			for (j = J + 1; j <= 2 * J; j++) {
				Wj = getAngel(tides.get(j - 1 - J));
				rslt[i][j] += Math.sin(Wj * ti) * Math.sin(Wi * ti);
			}
		}

	}

}
